function Pn = fn_u2Pn(n, H, T, un)

if abs(H/T)>2e-3

Pn = (1-exp(-un))/T ...
    .*2*((n+1)*coth(H/2/T*(n+1)) - coth(H/2/T)).*sinh(H/2/T)^2 ...
    ./(cosh(H/T*(n+1)) - exp(-un) - (1-exp(-un))*cosh(H/T));

else
    
    Pn = (1-exp(-un))*(n^2+2*n)/6/T*H/T./(n^2+2*n+exp(-un));
    
end

end